AD-A054  287 


UNCLASSIFIED 


AO 

*054287 


GENERAL  ELECTRIC  CO  SCHENECTADY  N V RESEARCH  AND  DEV— ETC  F/8  9/1 
TIME-DOMAIN  two-dimensional  COMruTER  SIMULATION  OF  THREE-TERMIN— ETC! 
1977  S P YU»  M TANTRAPORN  F44620-76-C-0043 

AFOSR-TR-78-0786  NL 


6 -78 


IDC  FILE  copy  AD  A054287 


ABOSE  TE-  7 8-  0 786 


-6  d c 


Time-Domain  Two-Dimensional  Computer  Simulation 
of  Three-Terminal  Semiconductor  Devices* 

S.P.  Yu+  and  W.  Tantraporn+ 


Abstract 


fTiI?l75[r1DD.DC 

MAY  28  1978 

15ED run 

B 


Special  schemes  for  treating  device  geometry  permit  the  use  of 
the  very  fast  Fourier  analysis  technique  for  solving  the  two  dimensional 
Poisson's  equation  so  that  the  dynamics  of  three  terminal  semiconductor 
devices  can  be  computer  simulated  economically.  FET's  and  bipolar  de- 
vices have  been  successfully  simulated. 

I.  Introduction 


In  recent  years,  computer  simulation  techniques  have  become  an  indis- 
pensable tool  in  solid-state  device  research  and  development.  This  is 
because  the  internal  physics  of  these  devices  and  their  interaction  with 
external  circuits  are  usually  governed  by  a system  of  highly  nonlinear 
partial  differential  equations.  The  equations  are  intractable  except  in 
certain  cases  when  they  can  be  linearized.  However,  the  linear  solution 
loses  some  important  aspects  of  the  device  physics  and  at  best  provides 
qualitative  information. 

One-dimensional  dynamic  computer  programs  provide  good  simulation  for 
many  devices  such  as  Gunn,  IMPATT  and  TRAPATT  diodes.  In  this  type  of 
program  the  algorithms  are  relatively  simple  to  formulate  and  the 
execution  time  is  relatively  short.  This  is  because  the  solution  of 
Poisson's  equation  can  be  reduced  to  the  inversion  of  a bidiagonal  matrix 
[2].  Three-terminal  devices  such  as  the  bipolar  transistor  and  the  FET, 
due  to  the  nature  of  their  physical  geometries  require  that  simulation 
must  account  for  at  least  two  dimensions.  The  algorithms  for  solving 
Poisson's  equation  in  two  dimensions  involves  inversion  of  a tridiagonal 
matrix  of  large  size.  The  usual  technique  for  computation  is  the  time 
consuming  over-relaxation  method.  Thus  for  reasonable  computer  process- 
ing time,  most  two  dimensional  semiconductor  device  simulations  are 
limited  to  the  static  case  [31  — [5]  .’ 


Through  judicious  choice  of  geometrical  representations  we  have  formu 
-lated  a dynamic  two-dimensional  program  for  the  simulation  of  three- 
terminal  semiconductor  devices.  For  bipolar  devices,  the  program  is  two- 
dimensional  in  the  collector  region  where  carrier  motion  transverse  to 
the  direction  of  the  applied  electric  field  has  a profound  effect  on  the 
device  physics.  The  emitter-base  junction,  owing  to  its  inherent 
geometry,  is  adequately  approximated  by  a one-dimensional  R-C  transmission 
line  with  a distributed  current  source.  This  scheme  for  treating  the 
device  geometry  permits  application  of  the  very  fast  Hockney's  technique 

[6]  to  two-dimensional  dynamic  simulation  of  the  transistor  and  the  CATT 

[7] ,  A description  of  the  bipolar  program  and  computational  results  are 
given  in  a companion  paper  [8].  Due  to  space  limitations,  this  paper  will 
focus  on  the  simulation  of  the  FET.  The  latter,  being  simpler  in  geometry, 
is  completely  simulated  dynamically  in  two  dimensions  as  described  below. 


* This  work  was  supported  by  the  US  Air  Force  Office  of  Scientific  Re- 
search, Contract  No.  F44620-76-C-0043 
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II.  Description  of  the  Computer  Program 

The  physical  models  for  FET  simulation  of  both  the  conventional  and 
the  double-sided  cases  are  shown  in  Fig.  1.  The  device  internal  physical 
effects  treated  in  our  program  are  the  following:  two-dimensional  carrier 
transport  resulting  from  electric  field  and  diffusion;  carrier  and  charge 
conservation  during  transport  and  avalanche  processes;  space  charge 
effects  resulting  from  mobile  and  fixed  charges;  field  and  time  dependent 
avalanche  process;  nonlinear  field-dependent  carrier  velocities,  and 
velocity  over-shoot  when  applicable  (See  Appendix) . The  external  bound- 
aries are:  Ohmic  contacts  on  the  source  and  drain  electrodes,  Schottky 
barrier  on  the  gate  electrodes,  electrode  voltages,  and  conservation  of 
terminal  currents. 

The  equation  governing  the  two  dimensional  device  are: 


Poisson's  equation 


and  carrier  transport 


V P 

?.  ( in  v ) + 
v " ' at 

v-  ip?) + If 
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The  symbols  used  in  the  above  equations  have  the  usual  meaning.  We 
proceed  now  to  qualitatively  discuss  the  procedure  for  solving  the  Poisson 
equation  (which  includes  manipulation  of  boundary  conditions),  the  con- 
tinuity equation  (which  includes  the  contact  problem  and  the  stability 
criterion)  and  the  definition  of  forces  on  carriers  and  their  velocities. 


The  device  is  partitioned  into  96  X 48  space  increments  (Ax,  Ay)  for 
the  solution  of  Poisson's  equation  and  carrier  transport,  as  shown  in  Fig. 
2.  For  computer  economy,  two  time  increments  At  and  Atfc  are  used  for  up- 
dating the  Poisson  equation  and  the  transport  eqSation  respectively.  The 
time  increments  are  related  by  an  integer  such  that  NAtfc  « At  . Thus  in 
the  duration  of  At  the  electric  distribution  remains  constan?  while  the 
carrier  motion  is  updated  N times  (N  = 2-*-12)  . 


The  use  of  Hockney's  technique  raises  some  complications  in  the  FET 
simulation  because  the  FET  contains  both  Dirichlet  and  Neumann  boundaries 
on  the  x » 0 and  48  edges,  while  the  Hockney  technique  is  applicable  to 
either  but  not  mixed  boundaries.  In  addition,  the  extent  of  the  boundary 
is  time  dependent.  In  order  to  alleviate  the  difficulties  the  Neumann 
boundaries,  which  arises  from  the  gate  depletion  region  (See  Fig.  2),  are 
converted  to  Dirichlet  in  each  time  increment.  This  could  be  done  because 
no  current  can  flow  across  the  exterior  boundaries  (x  = a)  of  the  deple- 
tion region  thus  the  condition  £$1  s q prevails.  This  condition,  together 

with  the  known  doping  density  and  potential  at  all  electrodes  permits  the 
evaluation  of  ^(a.y.t)  and  $2(a,y,t) . Thus  the  potential  <f>(a,y,t)  for 
o<y<b  is  completely  specified.  The  function  $(a,y,t)  is  Fourier  analyzed 
for  each  At  increment.  The  resultant  Fourier  components  are  then  applied 
to  approprlfte  location  in  Hockney's  algorithm.  The  vertical  boundaries 
at  y • o and  y - b are  considered  as  planes  of  even  symmetry.  This  is 
effected  by  using  only  a cosine  expansion  in  the  Fourier  analysis. 
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The  carrier  transport  Is  solved  by  an  explicit  technique.  Due  to  the 
negative  differential  mobility  in  GaAs,  formation  of  high  field  domain 
takes  place  under  certain  conditions.  Therefore,  it  is  necessary  to  use  a 
sufficiently  small  time  increment  size  Att  to  resolve  the  domain  formation 
and  decay.  The  condition  used  in  choosing  the  size  of  Att  to  ensure  a 
stable  solution  is  At.  v < Ur  . Further  refinements  in  solving  the  trans- 
A«  ' * 

port  equation  is  the  presentation  of  the  divergent  and  gradient  operators. 
Either  forward  or  backward  difference  is  used  depending  on  the  direction 
of  carrier  motion  and  diffusion  force.  This  technique  Improves  computa- 
tion accuracy  while  preserving  the  symmetry  (in  the  case  of  the  double- 
sided  FET) . The  ohmic  contact  at  the  source  electrode  is  ensured  by 
maintaining  the  mobile  carrier  density  in  the  interior  cells  adjacent  to 
the  electrodes  equal  to  the  background  doping.  The  surface  y m 0 may  be 
non-injecting  (Figure  la)  or  ohmic  (Figure  lb)  as  desired.  The  drain 
electrodes  and  the  surface  y = b act  as  an  infinite  sink  for  elections. 


Experimental  v-E  characteristic  for  GaAs  or  Si  is  used  in  specifying 
carrier  velocities  depending  on  the  device  material  being  simulated.  How- 
ever, instead  of  the  electric  field  alone  we  define  a force  which  includes 
diffusion: 


? * - 1 Vn 

f = -v4>  - Ivp 


for  electrons. 


for  holes. 


The  carrier  velocities  are: 


„v  = -v(|j|)  • b3/Lv| 
v 5 pv/|  V| 


The  forces  as  defined  are  consistent  with  the  concept  of  a quasi- 
Fermi  potential.  The  scalar  velocity  functions  are  those  available  in 
the  literature  as  v-E  characteristics.  For  electrons  in  GaAs  at  very  high 
frequencies  (>20  GHe) , the  effect  of  the  interband  delay  causing  the 
'velocity  overshoot' is  discussed  in  the  Appendix.  The  subroutine  dealing 
with  the  effect  has  been  checked  against  Shur's  calculation  [9],  but  in 
the  results  presented  below  the  subroutine  was  not  activated, 

III.  Results  and  Discussion 

Computational  results  indicate  that  the  96  X 48  cells  provide  good 
resolution  in  the  simulation  of  the  double-sided  device.  However,  they 
are  somewhat  inadequate  for  the  conventional  device.  This  is  because 
the  conventional  device  has  a semi- insulating  substrate  of  thickness  very 
large  compared  to  any  of  the  active  region  dimensions.  While  in  the  real  — 
world  the  potential  at  the  bottom  of  the  substrate  may  be  set  equal  to  _ * ' . 

zero  without  affecting  the  electric  field  distribution  in  the  device  Section 

active  region,  this  is  not  true  in  our  simulation  per  Figure  la.  The  “ • 

dimension  "a"  is  divided  into  48  equal  space  increments;  consequently  '•‘l0n  □ 

the  surface  x - 0 is  not  sufficiently  far  away  from  the  active  region  for  □ 

a good  approximation  as  an  equlpotential  surface.  The  same  problem  _ 

appears  to  exist  in  the  work  of  Barnes  and  Lomax  (10).  I 


1 Kimm/wMZ'jir  codes 
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The  speed  of  our  computer  program  is  approximately  two  orders  of  mag- 
nitude faster  than  those  using  the  relaxation  technique.  It  takes  about 
one  minute  to  complete  one  At  (which  consists  of  solving  the  Poisson's 
equation  once  the  transport  equation  N times)  in  the  96  X 48  mesh  points 
using  the  Honeywell-600  computer.  On  the  other  hand  it  takes  five  minutes 
to  compute  one  time  increment  in  a two  dimensional  1000  mesh-point  problem 
using  the  line  relaxation  technique  on  the  same  computer. 

Figure  4 shows  the  potential  distribution  in  a 'conventional'  FET, 
0.5y  thick  on  1 p insulator  substrate  grounded  at  the  lower  edge.  The 
electrode  arrangement  is  such  that  the  equipotential  contours  are  well 
displayed.  Because  of  the  poor  approximation  inherent  in  this  configura- 
tion discussed  above,  we  prefer  to  consider  the  'conventional'  FET  as  one 
half  of  the  double-sided  FET  and  simulate  only  the  latter. 

Figure  5 shows  the  potential  distribution  in  a double-sided  FET. 

Since  the  time,  length,  doping  density  follow  certain  scaling  laws  while 
computer  time  step  must  satisfy  the  stability  criterion,  a 2y  by  8y  space 
is  simulated.  Higher  frequency  operation,  except  the  velocity  overshoot 
effect,  can  be  scaled  from  the  result. 

Figure  6 shows  time-sequence  snapshots  of  the  electron  density  dis- 
tribution in  the  same  FET  at  various  times  indicated. 

Figure  7 shows  the  drain  current  as  a function  of  time  for  two 
values  of  the  gate  voltages.  In  the  Figures  5-*7  the  drain  and  gate 
voltages  are  applied  stiffly  at  t = 0.  Thus  the  response  time  of  the  FET 
is  being  studied.  The  peaks  and  valleys  in  the  drain  current  vs.  time 
suggest  the  Gunn  effect  in  the  FET,  and  this  is  evident  in  Figure  6 at 
times  40-*68  psec.  Note  that  there  are  more  than  100  data  points  for  each 
curve  in  Figure  7.  The  very  finely  sequenced  pictures  (not  included  here) 
strongly  suggest  two-dimensional  Gunn  domain  ripples  in  the  lower  field 
regions  in  the  source  and  drain,  which  might  manifest  themselves  as  very 
high  frequency  noises. 

In  conclusion,  we  have  qualitatively  discussed  our  computer  simula- 
tion for  two  dimensions,  and  have  presented  some  results  on  the  GaAs  FET. 
Large  signal  FET  performance  as  well  as  the  very  high  frequency  effect 
due  to  the  velocity  overshoot  will  be  carried  out.  Results  on  bipolar 
devices  are  reported  on  elsewhere. 
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APPENDIX 

An  approximate  treatment  of  the  'velocity  overshoot'  effect 

We  first  assume  that  if  there  was  only  the  light  mass  electronic 
energy  band  in  GaAs  the  velocity-field  v(E)  characteristic  would  be  such 
as  shown  in  Figure  8,  curve  a.  Curve  b in  the  Figure  represents  the 
equilibrium  v(E)  characteristic.  An  electron  having  a velocity  higher 
than  that  given  by  curve  b at  a constant  electric  field  E.  above  the 
threshold  field  E^may  be  assumed  to  decay  along  the  vertical  dashed 
line  shown  in  Figure  8 with  a time  constant  x(E^) 

v(t)  * V(E,)  + { V4(E.)  - V^E,))  e‘t/T(°  — l1) 


where  v^  is  for  curve  a,  v for  curve  b,  and  it  is  assumed  that  at  t * 0 
the  electron  velocity  is  given  by  Vj(E^). 


Consider  an  electric  field  as  a function  of  position  and  let  E » E ^ 
at  e = and  E>Et,  for  e>e  . An  electron  passing  i * « at  t * 0 going 
along  the  +e  direction  woul8  acquire  a velocity  v(b^)  at  b^  * bq  +4b. 


v(o  = vs  (£(*,)) ♦ [vf(F<i,>)  - vtwA]  **r  j 


i! 


l — d) 


And  at  e. 


b + IAb 
o 


v(*<)  * (*»>)  ** I H ( - vs( E C*4\)]  wpi- 

' «•!  1 


A? 


}-(.) 


The  denominator  in  the  exponential  term  should  contain  some  average 
values  of  v and  t between  the  two  adjacent  e positions.  The  choice  taken 
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in  (2)  and  (3)  is  to  facilitate  simple  updating  computation  of  v as  a 
function  of  a. 

Note  that  if  calculation  is  performed  along  the  a path  with  constant 
increment  Ae,  then  the  corresponding  At  will  inversely  vary  as  v(E(*)). 
This  is  in  fact  a convenient  way  to  calculate  the  effect  of  the  'velocity 
overshoot'  on  the  f of  the  transistor.  In  the  formulation  of  our  compu- 
ter simulation,  however,  the  space  grid  size  and  the  time  step  are  fixed, 
and  are  chosen  such  that  they  satisfy  the  stability  criterion  as  well  as 
minimizing  the  pseudo-diffusion  effect.  Thus  only  the  following  approxi- 
mate form  is  used  for  the  product  term  in  (3): 

where  K - A®/2  x 10?  cm/sec.  In  addition,  for  the  2 dimensional  simula- 
tion of  the  FET  a further  approximation  is  used  the  path  e is  along  the 
defined  y-direction  only. 

The  data  on  the  curve  a of  Figure  8 and  t(E)  are  obtained  from  the 
calculated  curves  published  by  Shur  (9) . It  is  found  that  the  curve  a 
can  be  represented  analytically  by: 

i.  s %0.43 

V jc  2.15  *10  + | E - 3*  io  i 

and  the  field  dependent  relaxation  time  is  given  by: 

-'3  . . -4  \ 1.2 

1 = 4 * 10  /i  E * 10  f 
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SOURCE  GATE  DRAIN 


Figure  5 Snapshots  of  potential  < 
double-sided  FET  with  the  gate  and 


SOURCE  GATE  DRAIN 


Figure  6 Time  sequence  (as  indicated,  in  psec)  snapshots  of  the  electron 
density  distribution  in  the  same  FET  as  in  Figure  5.  (The  grey  scale 
changes  between  the  "30"  and  "40"  frames).  Note  the  Gunn  oscillation, 
more  evident  as  the  shrinking  and  expanding  channel  width  near  the 
drain,  as  well  as  the  traveling  bulges  in  the  channel.  The  gate  voltage 
is  -4  volts. 
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Figure  7 The  drain  current  (passing  through  the  vertical  surface  beyond 
the  right  edge  of  the  gate  electrode)  as  a function  of  time,  for  two 
values  of  gate  voltages  on  the  same  FET  as  in  Figures  5 and  6. 


Figure  8 The  light-mass  velocity  branch  (curve  a)  and  the  equilibrium 
v(E)  characteristic  (curve  b)  of  electron  in  GaAs.  Curve  a as  well  as 
the  field  dependent  time  constant  of  decay  are  deduced  from  the  data 
published  by  Shur. 
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